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Among the ideas to be conveyed to students in an introductory quantum course, we have the 
pivotal idea championed by Dirac that functions correspond to column vectors (kets) and that 
differential operators correspond to matrices (ket-bras) acting on those vectors. The MATLAB 
(matrix-laboratory) programming environment is especially useful in conveying these concepts to 
students because it is geared towards the type of matrix manipulations useful in solving introductory 
quantum physics problems. In this article, we share MATLAB codes which have been developed at 
WPI, focusing on ID problems, to be used in conjunction with Griffiths' introductory text. 



Two key concepts underpinning quantum physics are 
the Schrodinger equation and the Born probability equa- 
tion. In 1930 Dirac introduced bra-ket notation for state 
vectors and operators.^ This notation emphasized and 
clarified the role of inner products and linear function 
spaces in these two equations and is fundamental to 
our modern understanding of quantum mechanics. The 
Schrodinger equation tells us how the state ^E" of a particle 
evolves in time. In bra-ket notation, it reads 



where H is the Hamiltonian operator and \^) is a ket or 
column vector representing the quantum state of the par- 
ticle. When a measurement of a physical quantity A is 
made on a particle initially in the state 5", the Born equa- 
tion provides a way to calculate the probability P{Ao) 
that a particular result Aq is obtained from the measure- 
ment. In bra-ket notation, it reads^ 

P{A,)^\{Aom? (2) 

where if \Ao) is the state vector corresponding to the 
particular result Ao having been measured, {Ao\ — |^o)^ 
is the corresponding bra or row vector and (^o|*) is thus 
the inner product between \Ao) and |^). In the Dirac 
formalism, the correspondence between the wavefunction 
vE'(f) and the ket I^E*) is set by the relation *(£) = {x\^), 
where \x) is the state vector corresponding to the particle 
being located at x. Thus we regard ^'(x) as a component 
of a state vector I^E"), just as we usuallj*^ regard ~ i-a 
as a component of a along the direction i. Similarly, we 
think of the Hamiltonian operator as a matrix 

H = Jd'j:\x)i^l-(^^^^\ix,t) + U{x)Yx\ (3) 

acting on the space of kets. 

While an expert will necessarily regard Eqs.(l-3) as a 
great simplification when thinking of the content of quan- 
tum physics, the novice often understandably reels under 
the weight of the immense abstraction. We learn much 
about student thinking from from the answers given by 
our best students. For example, we find a common error 



when studying ID quantum mechanics is a student treat- 
ing \E'(a;) and l^*) interchangeably, ignoring the fact that 
the first is a scalar but the ket corresponds to a column 
vector. For example, they may write incorrectly 



(p||*)|a;) = |*)(p|a;) (incorrect!) (4) 



or some similar abberation. To avoid these types of mis- 
conceptions, a number of educators and textbook authors 
have stressed incorporating a numerical calculation as- 
pect to quantum coursesi^i^iii^i^ The motive is simple. 
Anyone who has done numerical calculations can't help 
but regard a ket j^*) as a column vector, a bra (^| as a 
row vector and an operator H as a matrix because that is 
how they concretely represented in the computer. Intro- 
ducing a computational aspect to the course provides one 
further benefit: it gives the beginning quantum student 
the sense that he or she is being empowered to solve real 
problems that may not have simple, analytic solutions. 

With these motivations in mind, we have developed 
MATLAB codeslS' for solving typical 1 D problems found 
in the first part of a junior level quantum course based 
on Griffith's book^ We chose MATLAB for our pro- 
gramming environment because the MATLAB syntax is 
especially simple for the typical matrix operations used 
in ID quantum mechanics problems and because of the 
ease of plotting functions. While some MATLAB numeri- 
cal recipes have previously been published by others,— 
the exercises we share here are special because they em- 
phasize simplicity and quantum pedagogy, not numerical 
efficiency. Our point has been to provide exercises which 
show students how to numerically solve 1 D problems in 
such a way that emphasizes the column vector aspect of 
kets, the row vector aspect of bras and the matrix aspect 
of operators. Exercises using more efficient MATLAB 
ODE solvers or finite-element techniques are omitted be- 
cause they do not serve this immediate purpose. In part 
II of this article, we hope to share MATLAB codes which 
can be used in conjunction with teaching topics pertain- 
ing to angular momentum and non-commuting observ- 
ables. 
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I. FUNCTIONS AS VECTORS 

To start students thinking of functions as column 
vector-like objects, it is very useful to introduce them 
to plotting and integrating functions in the MATLAB 
environment. Interestingly enough, the plot command in 
MATLAB takes vectors as its basic input element. As 
shown in Program 1 below, to plot a function f{x) in 
MATLAB, we first generate two vectors: a vector of x 
values and a vector of y values where y — f{x). The com- 
mand plot{x,y/ r') then generates a plot window con- 
taining the points {xi.yi) displayed as red points (V). 
Having specified both x and y, to evaluate the definite 
integral ydx^ we need only sum all the y values and 
multiply by dx. 



%******************************************************* 

°/o Program 1: Numerical Integration cind Plotting using MATLAB 

N=1000000; % No. of points 

L=500; °/o Range of x: from -L to L 

x=linspaceC-L,L,K) ' ; % Generate column vector with N 

% X values ranging from -L to L 
dx=xC2)-x(l) ; °/o Distance between adjacent points 

7o Alternative Trial functions: 

°/a To select one, take out the comment commEind % at the beginning. 
°/oy=expC-x. "2/16) ; % Gaussicin centered at x=0 

yoy=(C2/pi)~0.5)*exp(-2*(x-l) . "2) ; 7. Normed Gaussian at x=l 
yoy=(l-x . ~2) . ~-l ; % Symmetric fen which blows up at x=l 

yoy=(cos (pi*x) ) . "2 ; °/o Cosine fen 

yoy=expCi*pi*x) ; % Complex exponential 

yy=sinCpi*x/L) . *cos (pi*x/L) ; y. Product of sinx times cosx 
yy=sinCpi*x/L) . *sin(pi*x/L) ; y. Product of sinCnx) times sin(mx) 
yA=100; y=sinCx*A) ./(pi*x) ; y. Rep. of delta fen 

A=20; y=(sin(x*A) .'2) ./Cpi*CA*x.'2)) Another rep. of delta fen 
y Observe: numerically a function y(x) is represented by a vector! 
y Plot a vector/function 

plotCx,y); y Plots vector y vs . x 

yplot(x,realCy) , 'r ' , x, imag(y) , 'b'); "L Plots real&imag y vs. x 
axis C [-2 2 7]); y Optimized axis parameters for sinx~2/pix~2 
yaxis([-2 2 -8 40]); "L Optimized axis parameters for sinx/pix 

y Numerical Integration 

sum(y)*dx 'L Simple numerical integral of y 

trapz(y)*dx % Integration using trapezoidal rule 



II. DIFFERENTIAL OPERATORS AS 
MATRICES 

Just as f{x) is represented by a column vector |/) in 
the computer, for numerical purposes a differential opera- 
tor D acting on f{x) is reresented by a matrix D that acts 
on I/). As illustrated in Program 2, MATLAB provides 
many useful, intuitive, well-documented commands for 
generating matrices D that correspond to a given I?.— 
Two examples are the commands ones and diag. The 
command ones{a,b) generates an a x 6 matrix of ones. 
The command diag{A, n) generates a matrix with the el- 
ements of the vector A placed along the n^^ diagonal and 
zeros everywhere else. The central diagonal corresponds 
to n = 0, the diagonal above the center one corresponds 



to n = 1, etc.). 

An exercise we suggest is for students to verify that 
the derivative matrix is not Hermitian while the deriva- 
tive matrix times the imaginary number i is. This can 
be very valuable for promoting student understanding if 
done in conjunction with the proof usually given for the 
differential operator. 



7o Program 2: Calculate first and second derivative numerically 
°/o showing how to write differential operator as a matrix 

7o Parameters for solving problem in the interval < x < L 

L = 2*pi ; °/o Interval Length 

N = 100; °/o Ko . of coordinate points 

X = linspace (0 , L , K) ' ; % Coordinate vector 

dx = x(2) - x(l); °/o Coordinate step 

°/„ Two-point finite-difference representation of Derivative 
D=(diag(ones((N-l) ,1) ,l)-diag(ones((N-l) , 1) , -1) )/ (2*dx) ; 
y. Next modify D so that it is consistent with f(0) = fCL) = 
D(l,l) = 0; D(l,2) = 0; D(2,l) = 0; 7. So that f(0) = 
D(K,K-1) = 0; D(K-1,K) = 0; D(K,K) = 0; 7. So that f(L) = 

7 Three-point finite-difference representation of Laplacian 
Lap = (-2*diag(ones(N,l) ,0) + diagCones ( (N-1) , 1) , 1) ... 

+ diagCones ((K-1) ,1) ,-l))/(dx-2) ; 
7 Next modify Lap so that it is consistent with f(0) = f (L) = 
Lap(l,l) = 0; Lap(l,2) = 0; Lap(2,l) = 0; '/. So that f(0) = 
Lap(N,K-l) = 0; Lap(N-l,N) = 0; LapCN.N) = 0;7 So that f(L) = 

7 To verify that D*f corresponds to taking the derivative of f 
7 and Lap*f corresponds to taking a second derviative of f, 
7 define f = sinCx) or choose your own f 
f = sinCx) ; 

7 And try the following: 

Df = D*f; Lapf = Lap*f; 

plotCx,f, 'b' ,x,Df, 'r' , x,Lapf, 'g') ; 

axisC[0 5 -1.1 1.1]); 7 Optimized axis parameters 

7 To display the matrix D on screen, simply type D and return . . . 
D 7 Displays the matrix D in the workspace 
Lap 7 Displays the matrix Lap 



III. INFINITE SQUARE WELL 

When solving Eq. ([1]) , the method of separation of vari- 
ables entails that as an intermediate step we look for the 
separable solutions 

\^E{t)) = \^Emexp{-iEt/h) (5) 

where |^£;(0)) satisfies the time-independent Schrodinger 
equation 

i/|*i5(o)) = |*b(o)). (6) 

In solving Eq. ^ we are solving for the eigenvalues E 
and eigenvectors |^'_e(0)) of H. In MATLAB, the com- 
mand [V, E] = eig{H) does precisely this: it generates 
two matrices. The first matrix V has as its columns the 
eigenvectors \^e{^))- The second matrix _E is a diago- 
nal matrix with the eigenvalues Ei corresponding to the 
eigenvectors \^ eA^)) placed along the central diagonal. 
We can use the command E = diag{E) to convert this 
matrix into a column vector. In Program 3, we solve for 
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% Program 4: Find several lowest eigenmod.es V(x) and 
7o eigenenergies E of ID Schrodinger equation 
t -l/2*hbar-2/mCd2/dx2)V(x) + UCx)VCx) = EVCx) 

7o for arbitrary potentials UCx) 

°/o Parameters for solving problem in the interval -L < x < L 
% PARAMETERS: 

L = 5; °/o Interval Length 

N = 1000; % Ko of points 

X = linspace(-L,L,N) ' ; % Coordinate vector 

dx = x(2) - x(l); °/o Coordinate step 

% POTENTIAL , choose one or make your own 

U = l/2*100*x . " (2) ; 7o quadratic harmonic oscillator potential 
°/,U = l/2+x.~(4); 7. quartic potential 

7o Finite square well of width 2w cind depth given 
Zw = L/50; 

7U = -500* (heaviside (x+w) -heaviside Cx-w) ) ; 

7 Two finite square wells of width 2w and distance 2a apart 
7w = L/50; a=3*w; 

7U = -200* (heaviside (x+w-a) - heaviside (x-w-a) ... 

7 + heaviside Cx+w+a) - heaviside Cx-w+a) ) ; 

7 Three-point finite-difference representation of Laplacian 
7 using sparse matrices, where you save memory by only 
7 storing non-zero matrix elements 

e = ones(K,l); Lap = spdiags([e -2*e e] , [-1 1] ,N,N)/dx~2; 
7 Total Hamiltonian 

hbar = 1; m = 1; 7 constants for Hamiltonian 

H = -l/2*Chbar"2/m)*Lap + spdiagsCU,0,N,N) ; 

7 Find lowest nmodes eigenvectors and eigenvalues of sparse matrix 
nmodes = 3; options. disp = 0; 

[V,E] = eigs (H, nmodes sa' , options) ; 7 find eigs 

[E,ind] = sort (diagCE) ) ;7 convert E to vector and sort low to high 

V = VC:,ind); 7 rearrange corresponding eigenvectors 

7 Generate plot of lowest energy eigenvectors VCx) and U(x) 

Use = U*maxCabs CVC : ) ) ) /max(abs (U) ) ; 7 rescale U for plotting 

plot Cx,V,x, Use , ' — k'); 7 plot V(x) and rescaled UCx) 

7 Add legend showing Energy of plotted V(x) 

lgnd_str = [repmatC'E = ' ,nmodes , 1) ,num2str CE)] ; 

legendClgnd_str) 7 place lengend string on plot 

shg 



the eigenfunctions and eigenvalues for the infinite square 
well Hamiltonian. For brevity, we omit the commands 
setting the parameters L, A^, and dx. 



7******************************************************* 

7 Program 3: Matrix representation of differential operators, 
7 Solving for Eigenvectors & Eigenvalues of Infinite Square Well 

7 For brevity we omit the commainds setting the parameters L, N, 
7 X cind dx; We also omit the commainds defining the matrix Lap. 
7 These would be the same as in Program 2 above. 

7 Total Hamiltonian where hbar=l cind m=l 
hbar = 1; m = 1; H = - (1/2) * (hbar-2/m) *Lap ; 

7 Solve for eigenvector matrix V cind eigenvalue matrix E of H 
[V,E] = eigCH) ; 

7 Plot lowest 3 eigenfunctions 

plotCx,VC: ,3) , 'r' ,x,V(: ,4) , 'b' ,x,VC: ,5) , 'kO ; shg; 
E 7 display eigenvalue matrix 

diagCE) 7 display a vector containing the eigenvalues 



Note that in the MATLAB syntax the object y(:,3) 
specifies the column vector consisting of all the elements 
in column 3 of matrix V . Similarly 1^(2, :) is the row 
vector consisting of all elements in row 2 of V\ y(3, 1) is 
the element at row 3, column 1 of and F(2, 1 : 3) is 
a row vector consisting of elements F(2, 1), 1/(2,2) and 

y(2,3). 

IV. ARBITRARY POTENTIALS 

Numerical solution of Eq. ([T]) is not limited to any par- 
ticular potential. Program 4 gives example MATLAB 
codes solving the time independent Schrodinger equa- 
tion for finite square well potentials, the harmonic os- 
cillator potential and even for potentials that can only 
solved numerically such as the quartic potential U = x^. 
In order to minimize the amount of RAM required, the 
codes shown make use of sparse matrices, where only 
the non-zero elements of the matrices are stored. The 
commands for sparse matrices are very similar to those 
for non-sparse matrices. For example, the command 
[V, -E] = eigs{H, nmodes..) provides the nmodes lowest 
energy eigenvectors V of of the sparse matrix H. 

Fig. [T] shows the plot obtained from Program 4 for the 
potential U = ^ - 100 • x^. Note that the 3 lowest energies 
displayed in the figure are just as expected due to the 
analytic formula 

E = nu;(^n+^^ (7) 
with n — integer and a; — -\/— = 10 rad/s. 



V. A NOTE ON UNITS IN OUR PROGRAMS 

When doing numerical calculations, it is important to 
minimize the effect of rounding errors by choosing units 
such that the variables used in the simulation are of the 



order of unity. In the programs presented here, our fo- 
cus being undergraduate physics students, we wanted to 
avoid unnecessarily complicating matters. To make the 
equations more familiar to the students, we explicitly left 
constants such as Yi in the formulas and chose units such 
that ti = \ and to = 1. We recognize that others may 
have other opinions on how to address this issue. An al- 
ternative approach used in research is to recast the equa- 
tions in terms of dimensionless variables, for example by 
rescaling the energy to make it dimensionless by express- 
ing it in terms of some characteristic energy in the prob- 
lem. In a more advanced course for graduate students or 
in a course in numerical methods, such is an approach 
which would be preferable. 



VI. TIME DEPENDENT PHENOMENA 

The separable solutions \^Eify) are only a subset of 
all possible solutions of Eq. ((!]). Fortunately, they are 
complete set so that we can construct the general solution 
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FIG. 1: Output of Program 4, which plots the energy eigen- 
functions V{x) and a scaled version of the potential U{x) — 
1/2 • 100 • x'^ . The corresponding energies displayed within 
the figure legend, 4.99969, 14.9984 and 24.9959, are, within 
rounding error, precisely those expected from Eq. ([7]) for the 
three lowest-energy modes. 



via the linear superposition 

1^(0) - aE\-9EiO))expi~iEt/h) 



(8) 



where are constants and the sum is over all possible 
values of E. The important difference between the sep- 
arable solutions jl]) and the general solution ([8|) is that 
the probability densities derived from the general solu- 
tions are time-dependent whereas those derived from the 
separable solutions are not. A very apt demonstration 
of this is provided in the Program 5 which calculates the 
time-dependent probability density p{x, t) for a particle 
trapped in a a pair of finite-square wells whose initial 
state |5'(0)) is set equal to the the equally- weighted su- 
perposition 



|vI/(0)) = -i=(|vI/s„) + |*S,)) 



(9) 



of the ground state I^'bq) and first excited state \^ Ei) 
of the double well system. As snapshots of the program 
output show in Fig. [21 the particle is initially completely 
localized in the rightmost well. However, due to Eq ^ E\, 
the probability density 



p(x,t) = i[|^B„(x)|2 + \^EA^)f + 

2\^e,,(xWeA^)\cos' {{El - E2)t/h) ; 



is time-dependent, oscillating between the p{x) that cor- 
responds to the particle being entirely in the right well 

p{x)^\\4>Eo{x)\ + \4^eAx)\\^ (11) 
and p{x) for the particle being entirely in the left well 



p{x) = I \ikEo{x)\ - \i^EAx)\ I ■ 



(12) 



By observing the period with which p(x, t) oscillates 
in the simulation output shown in Fig. [2] students can 
verify that it is the same as the period of oscillation 
2TTh/{Ei - E2) expected from Eq. (|T0| . 



°/o Program 5: Calculate Probability Density as a function of time 
°/o for a particle trapped in a double-well potential 

7o Potential due to two square wells of width 2w 
7o and a distance 2a apart 
w = L/50; a = 3*w; 

U = -100+C heaviside Cx+w-a) - heaviside (x-w-a) ... 

+ heaviside (x+w+a) - heaviside (x-w+a) ) ; 

°l Finite-difference representation of Laplacian and Hamiltonian, 
7o where hbar = m = 1 . 

e = ones(K,l); Lap = spdiagsCEe -2*e e] , [-1 1] ,N,N)/dx~2; 
H = -Cl/2)*Lap + spdiagsCU,0,N,N) ; 

7 Find and sort lowest nmodes eigenvectors and eigenvalues of H 
nmodes = 2; options. disp = 0; [V,E] = eigs (H , nmodes sa' , options) ; 
[E,ind] = sort (diagCE) ) convert E to vector and sort low to high 
V = VC:,ind); % rearrange coresponding eigenvectors 

7, Rescale eigenvectors so that they are always 
7o positive at the center of the right well 
for c = 1: nmodes 

V(:,c) = VC: ,c)/sign(VCC3*N/4) ,c)) ; 

end 

7 Compute and display normalized prob. density rhoCx,T) 

7 Parameters for solving the problem in the interval < T < TP 
TP = 10; 7 Length of time interval 

NT = 100; 7 No. of time points 

T = linspace(0,TF,NT) ; 7 Time vector 

7 Compute probability normalisation constant (at T=0) 
psi_o = 0.5*VC : ,l)+0.5+V( : ,2) ; 7 wavefunction at T=0 
sq_norm = psi_o ' *psi_o*dx ; 7 square norm = |<ff |ff>|~2 



Use = U*maxCabsCVC:)))/max(abs(U)) ; 



7 rescale U for plotting 



7 Compute cind display rhoCx,T) for each time T 

for t=l:NT; 7 time index parameter for stepping through loop 
7 Compute wavefunction psi(x,T) and rhoCx,T) at T=TCt) 
psi = 0.5*VC: ,l)*exp(-i*ECl)*T(t)) ... 

+ 0.5*V(: ,2)*expC-i*E(2)*TCt)); 
rho = conj (psi) . *psi/sq_norm; 7 normalized probability density 

7 Plot rho(x,T) and rescaled potential energy Use 

plot (x, rho, 'o-k' , x, Usc,'.-b'); axis([-L/8 L/8 -1 6]); 

lgnd_str = [repmatCT = ' , 1 , 1) ,num2str (T) ] ; 

text (-0.12,5.5, lgnd_str , ' FontSize ' , 18) ; 

xlabelCx [m] ' , 'PontSize', 24); 

ylabel ( 'probability density [1/m] ' , ' FontSize ' , 24) ; 



pause (0 . 05) ; 
end 



7 wait 0.05 seconds 



(10) 
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lates and displays the time evolution of a wavepacket for 
one of two possible potentials, either J7 = or a step 
potential U = UoQ{x — L). The initial wavepacket is 
generated as the Fast Fourier Transform of a Gaussian 
momentum distribution centered on a particular value 
of the wavevector fco- Because the wavepacket is com- 
posed of a distribution of different ks, different parts of 
the wavepacket move with different speeds, which leads 
to the wave packet spreading out in space as it moves. 

While there is a distribution of velocities within the 
wavepacket, two velocities in particular are useful in char- 
acterizing it. The phase velocity = oj/k = E/p = 
hko/ 2-171 is the velocity of the plane wave component 
which has wavevector ko- The group velocity Vg = hko/m 
is the velocity with which the expectation value < a; > 
moves and is the same as the classical particle velocity as- 
sociated with the momentum p = hk. Choosing U ~ 0, 
students can modify this program to plot < x > vs t. 
They can extract the group velocity from their numerical 
simulation and observe that indeed Vg = 2vw for a typ- 
ical wave packet. Students can also observe that, while 
Vg matches the particle speed from classical mechanics, 
the wavepacket spreads out as time elapses. 

In Program 6, we propagate the wave function forward 
via the formal solution 

\^{t))=exp{~im/h)\^{0)), (13) 




-0.6-0.4-0.2 0.2 0.4 0.6 
X [m] 



FIG. 2: The probability density p(x,t) output from Program 
5 for a particle trapped in a pair of finite-height square well 
potentials that are closely adjacent. The initial state of the 
particle is chosen to be |^'(0)) ~ |i5o) -I- \E\). Shown is p{x,t) 
plotted for times T = {0, 0.25, 0.5, 0.75, 1.0} x 2-Kh/{Ei - E2). 



VII. WAVEPACKETS AND STEP POTENTIALS 

Wavepackets are another time-dependent phenomenon 
encountered in undergraduate quantum mechanics for 
which numerical solution techniques have been typically 
advocated in the hopes of promoting intuitive acceptance 
and understanding of approximations necessarily invoked 
in more formal, analytic treatments. Program 6 calcu- 



where the Hamiltonian matrix H is in the exponential. 
This solution is equivalent to Eq. ([5]) , as as can be shown 
by simple substitution. Moreover, MATLAB has no trou- 
ble exponentiating the matrix that numerically repre- 
senting the Hamiltonian operator as long as the matrix 
is small enough to fit in the available computer memory. 

Even more interestingly, students can use this method 
to investigate scattering of wavepackets from various po- 
tentials, including the step potential U = UoQ{x — L/2). 
In Fig. [21 we show the results of what happens as the 
wavepacket impinges on the potential barrier. The pa- 
rameters characterizing the initial wavepacket have been 
deliberately chosen so that the wings do not fall outside 
the simulation area and initially also do not overlap the 
barrier on the right. If (E) <C Uo, the wavepacket is 
completely reflected from the barrier. If (E) ^ Uo, a 
portion of the wave is is reflected and a portion is trans- 
mitted through. If (E) ^ Uo, almost all of the wave is 
transmitted. 

In Fig. [5] we compare the reflection probability R cal- 
culated numerically using Program 6 with R calculated 
by averaging the single-mode^"'^ expression 



RiE) = 



VE-^E-Uo 



y/E + VE-Uo 



(14) 



over the distribution of energies in the initial wavepacket. 
While the numerically and analytically estimated R are 
found to agree for large and small {E)/Uo, there is a 
noticeable discrepancy due to the shortcomings of the 
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°/o Program 6: Wavepacket propagation using exponential of H 

Parameters for solving the problem in the interval < x < L 
L = 100; ■/. Interval Length 

N = 400; •/. No of points 

X = linspace (0 , L , K) ' ; X Coordinate vector 

dx = x(2) - x(l); °/o Coordinate step 

% Parameters for making intial momentum space wavef unction phi(k) 

ko = 2; % Peak momentum 

a = 20; % Momentum width parameter 

dk = 2*pi/L; % Momentum step 

km=N*dk; % Momentum limit 

k=linspaceCO,+km,N) ' ; % Momentum vector 

°/o Make psi(x,0) from Gaussian kspace wavefunction phi(k) using 
7o fast fourier transform : 

phi = expC-a* Ck-ko) . '2) . *exp(-i*6*k . '2) ; 'h unnormalized phi Ck) 
psi = ifft(phi); % multiplies phi by expikx and integrates vs. x 
psi = psi/sqrt (psi ' *psi*dx) ; % normalize the psiCx,0) 

7o Expectation value of energy; e.g. for the parameters 
7. chosen above <E> = 2.062. 

avgE = phi ' *0 . 5*diag(k. "2,0)*phi*dk/(phi' *phi*dk) ; 
7 CHOOSE POTENTIAL UCX) : Either U = OR 

7 U = step potential of height avgE that is located at x=L/2 

7U = 0*heaviside (x- CL/2) ) ; 7 free particle wave packet evolution 

U = avgE*heaviside (x- CL/2) ) ; 7 scattering off step potential 

7 Finite-difference representation of Laplacian and Hamiltonian 
e = ones(N,l); Lap = spdiags([e -2*e e] , [-1 1] ,N,N)/dx"2; 
H = -(l/2)*Lap + spdiags(U,0,N,N) ; 

7 Parameters for computing psi(x,T) at different times < T < TP 

NT = 200; 7 No. of time steps 

TF = 29; T = linspace(0,TF,NT) ; 7 Time vector 

dT = T(2)-T(l); 7 Time step 

hbar = 1; 

7 Time displacement operator E=expC-iHdT/hbar) 

E = expm(-i*f ull (H) *dT/hbar) ; 7 time desplacement operator 

7 Simulate rho(x,T) and plot for each T 

for t = 1:NT; 7 time index for loop 

y, calculate probability density rhoCx,T) 
psi = E*psi; 7 calculate new psi from old psi 

rho = conj (psi) . +psi ; 7 rho(x,T) 

plot (x ,rho , 'k' ) ; 7 plot rho(x,T) vs. x 

axis([0 L 0.15]); 7 set x,y axis parameters for plotting 
xlabelCx [m] ' , 'FontSize', 24); 

ylabel ( 'probability density [1/m] FontSize ' , 24); 



pause CO . 05) ; 



7 pause between each frame displayed 



end 



7 Calculate Reflection probability 
R=0; 

for a=l:N/2; 

R=R+rho(a) ; 

end 
R=R*dx 



numerical simulation for {E)/Uo ~ 1. This discrepancy 
can be reduced significantly by increasing the number of 
points in the simulation to 1600 but only at the cost of 
significantly slowing down the speed of the computation. 
For our purposes, the importance comparing the analyt- 
ical and numerical calculations is that it gives student 
a baseline from which to form an opinion or intuition 
regarding the accuracy of Eq. . 




10 20 30 40 50 60 70 80 90 100 
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FIG. 3: Output of Program 6 showing a wavepacket encoun- 
tering step potential of height ~ (E) located a,t x/L = 0.5 at 
different times. 



VIII. CONCLUSIONS 

One benefit of incorporating numerical simulation into 
the teaching of quantum mechanics, as we have men- 
tioned, is the development of student intuition. Another 
is showing students that non-ideal, real-world problems 
can be solved using the concepts they learn in the class- 
room. However, our experimentation incorporating these 
simulations in quantum physics at WPl during the past 
year has shown us that the most important benefit is a 
type of side-effect to doing numerical simulation: the ac- 
ceptance on an intuitive level by the student that func- 
tions are like vectors and differential operators are like 
matrices. While in the present paper, we have only had 
sufficient space to present the most illustrative MATLAB 
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codes, our goal is to eventualy make available a more 
complete set of polished codes is available for download- 
ing either from the authors or directly from the file ex- 
change at MATLAB Central^ 



FIG. 4: The reflection probability R vs. {E)/Uo- The dashed 
line is simply Eq. (|14|l where we subsitute E = {E), the solid 
line is Eq. (|14|) averaged over the energy distribution in the 
incident wavepacket, and the points are numerical results ob- 
tained using Program 6, where the horizontal distance be- 
tween points is oeIUo where as is the standard deviation of 
the energy distribution in the wavepacket. 
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